#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "CGSolver.h"
#include "clapack.h"

int i,jmt;

void cgmethod(double A[], double b[], double x[],int iz,int lda){
    integer M = iz;
    integer N = iz;
    integer LDA = lda;
    address  T = 'n';
    double Ap[iz];
    double pAp;
    double alpha;
    double rsnew;
    double rold[iz];
    integer Kl = 2;
    integer Ku = 1;
    double al = -1;
    double be = 1;
    integer incx = 1;
    integer incy = 1;
    double zero = 0;
    //b = b-A*x0
    //r = b

    dgbmv_(&T,&M,&N,&Ku,&Kl,&al,A,&LDA,x,&incx,&be,b,&incy);
    //p = b
    double p[iz];
    dcopy_(&M,b,&incx,p,&incy);
    //rsoldr'*r
    double rsold;
    rsold = ddot_(&M,b,&incx,b,&incy);
    for(i = 0;i<10000;i++){
        //Ap = A * p


        dgbmv_(&T,&M,&N,&Ku,&Kl,&be,A,&LDA,p,&incx,&zero,Ap,&incy);
       /* for(jmt=0;jmt<4;jmt++){
            printf("value AP %f \n",Ap[jmt]);
        }*/
        pAp = ddot_(&M,p,&incx,Ap,&incy);
        //printf("value pAp %f \n",pAp);
        alpha = rsold/pAp;
        //printf("value alpha %f \n",alpha);
        daxpy_(&M,&alpha,p,&incx,x,&incy);
        /*for(jmt=0;jmt<4;jmt++){
            printf("value x %f \n",x[jmt]);
        }*/
        double min_alp = -alpha;
        daxpy_(&M,&min_alp,Ap,&incx,b,&incy);
        /*for(jmt=0;jmt<4;jmt++){
            printf("value new b %f \n",b[jmt]);
        }*/
        rsnew = ddot_(&M,b,&incx,b,&incy);

        if (sqrt(rsnew)<pow(10,-10)){
        printf("iter CG %d \n",i);
            break;
        }
        //memcpy(rold,b,m*sizeof(double));
        dcopy_(&M,b,&incx,rold,&incy);
        double rs = rsnew/rsold;
        daxpy_(&M,&rs,p,&incx,b,&incy);
        //memcpy(p,b,m*sizeof(double));
        dcopy_(&M,b,&incx,p,&incy);
        /*for(jmt=0;jmt<4;jmt++){
            printf("value new p %f \n",p[jmt]);
        }*/
        //memcpy(b,rold,m*sizeof(double));
        dcopy_(&M,rold,&incx,b,&incy);
        /*for(jmt=0;jmt<4;jmt++){
            printf("value new b %f \n",b[jmt]);
        }*/
        rsold = rsnew;



    }






}
